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The neural dynamics of the nematode C. elegans are experimentally low-dimensional and cor¬ 
respond to discrete behavioral states, where previous modeling work has found neural proxies for 
some of these states. Experimental results further suggest that dynamics may be understood as long- 
timescale transitions between multiple low-dimensional attractors. To identify multistable regimes 
of our model, we develop a method for systematic generation of bifurcation diagrams and their 
analysis in an interpretable low-dimensional subspace, showing the existence and nature of multi¬ 
stable input responses at a glance. Stimulation of the PLM neuron pair, experimentally associated 
with forward movement and shown in simulation to drive a limit cycle, defines our low-dimensional 
projection space. We then obtain bifurcation diagrams for single-neuron excitation over a range of 
amplitudes and which classify whether the dynamics in this projection space are associated with 
a limit cycle, fixed point, or multiple states. In the specific case of compound input into both 
the PLM pair and ASK pair we discover bistability of a limit cycle and a fixed point, with tran¬ 
sitional timescales between different states being much longer than other timescales in the system. 

This suggests consistency of our model with the characterization of dynamics in neural systems as 
long-timescale transitions between discrete, low-dimensional attractors corresponding to behavioral 
states. Our methodology thus prescribes a method for identifying these states and transitions in 
response to arbitrary input. 


I. INTRODUCTION 

Understanding the functional responses and control of 
high-dimensional networked dynamical systems is of crit¬ 
ical importance across the physical, engineering and bi¬ 
ological sciences. In many such systems, even with large 
numbers of nonlinearly interacting nodes, meaningful in¬ 
put and output are dominated by low dimensional spatio- 
temporal patterns of activity m- Indeed, the underly¬ 
ing networked dynamics can be thought to be dominated 
by trajectories that evolve on low-dimensional attractors 
and/or induced transient trajectories between attractors. 
As a specihc biophysical example, neuronal networks, 
which are typified by high-dimensional networks of neu¬ 
rons, display robust functional responses and behavioral 
assays that are encoded by such low-dimensional attrac¬ 
tors or transient trajectories HHII]. 

The nematode Caenorhahditis elegans is an impor¬ 
tant model system in understanding how these neu¬ 
ronal networks generate robust functional responses to 
inputs, partly due to the fact that the connectivity be¬ 
tween its 302 neurons (its connectome) has been re¬ 
solved [nna. C. elegans is capable of a wide range of 
behaviors over various timescales na, yet experimental 
studies suggest that these behaviors are fundamentally 
low-dimensional m, and the behaviors of the worm can 
be understood as low-dimensional trajectories on attrac¬ 
tors between which it will transition stochastically [Ill- 

While the exact role of the connectome in neuronal 
computation remains unresolved and in general contro¬ 
versial, it has been shown that simple models of C. el¬ 
egans neural dynamics (combining specific connectivity 
data with simple unfit parameter estimates and dynam¬ 


ics) are capable of generating non-trivial, qualitatively 
correct responses to given stimuli m- This suggests 
that such computational modeling can be informative in 
understanding how the system generates behavioral re¬ 
sponses. It is therefore of interest to consider whether or 
not models that are capable of producing neural proxies 
for behavioral responses (as in m) are further capable 
of characterizing experimentally observed multistable at¬ 
tractor dynamics. 

One motive in the search for multistability is that the 
previous study in m finds a neural proxy for behavior 
consisting simply of a single limit cycle within the system. 
On the other hand, experimental evidence suggests that 
many neural responses are better described as transient 
trajectories between multiple attractors [7], rather than 
the traditional dynamical systems view in which behavior 
is described by dynamics on a stable attractor. Within 
this paradigm it is important not only that multistabil¬ 
ity should exist, but that transients with biophysically 
relevant behavioral timescales should exist. 

In this paper, we explore the input space of our full 
model for the neuronal network dynamics of C. elegans, 
developed in m, and find that various multistabilities 
arise in response to inputs. Performing direct neuronal 
simulations to reveal such transitions is a formidable 
task, since the input space is large and neuronal simu¬ 
lations produce high-dimensional outputs which are dif¬ 
ficult to interpret. We therefore develop a systematic 
methodology to explore responses to complex inputs and 
understand the dynamics within a framework of low¬ 
dimensional attractor dynamics. Specifically, for a cho¬ 
sen input vector we generate a bifurcation diagram (using 
the amplitude of the constant-in-time input as our bifur- 
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cation parameter) to identify multistability. We use an 
interpretable low-dimensional projection (as defined by 
forward motion) to characterize the dynamics between 
multiple attractors as identified by the bifurcation dia¬ 
grams. 

With this framework we survey all input vectors cor¬ 
responding to single-neuron current injections and find 
bifurcations corresponding to limit cycles and multiple 
attractors. For some of these input vectors which induce 
multistability, simulated transient dynamics are on much 
slower (on the order of seconds to tens of seconds) than 
any intrinsic neuronal timescales (which in our model do 
not exceed a few hundred milliseconds). The transient 
trajectories themselves are low-dimensional and could be 
associated with network-produced functionalities, such as 
neural proxies for movement. Thus our simulations of 
connectomic dynamics are in agreement with behavioral 
observations of C. elegans and help support recent bio¬ 
physical conjectures that the transients themselves are 
critical in understanding behavioral assays [7]. 

As a particular example, we choose input into the PLM 
neuron pair, which is known experimentally to excite for¬ 
ward motion m and within our model creates a two- 
dimensional limit cycle response m- We then use the 
low-dimensional PLM response plane to consider the dy¬ 
namics of a compound input vector PLM-I-ASK, where 
ASK stimulation is known to facilitate transitions (i.e. 
turns) [18) . Our bifurcation analysis reveals that this in¬ 
duces bi-stability, in which the system goes either into 
a fixed point or a limit cycle. Transient timescales are 
shown to be considerably longer in this bistable case than 
the intrinsic timescales of the system. This allows for 
long timescales in the system arising from transitions be¬ 
tween discrete, low-dimensional attractors corresponding 
to behavioral states, consistent with the experimentally- 
based framework m- This input scenario demonstrates 
how our bifurcation analysis methodology prescribes a 
generic approach for identifying multi-stable states and 
transitions between them in response to arbitrary inputs. 
Since we model neurons as identical except for their con¬ 
nectivity, it further indicates that their connectivity alone 
can encode the creation and destruction of multiple be¬ 
havioral attractors and transitions between them. 


II. SIMULATION OF C. ELEGANS 
CONNECTOME DYNAMICS 

A. Model for Coupled Neural Dynamics 

The dynamic model used here is constructed to repre¬ 
sent the graded responses of the neurons of C. elegans. 
Experiments show that many neurons in the organism 
are effectively isopotential, so that the membrane volt¬ 
age can be used as a state variable for characterizing 
the neuron HD. Wicks et al. [20] used this to construct 
a single-compartment membrane model for neuron dy¬ 
namics. Building on these findings in El, we were able 
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FIG. 1. Voltage dynamics of forward-motion motorneurons 
(neurons of classes DB, DD, VB and VD) in response to the 
following sensory inputs; in Panel A, an input of 2 x 10"^ 
(Arb. Units) into the PLM sensory neuron pair (known ex¬ 
perimentally to drive forward motion El) ; in Panel B, an 
input of 2 X 10"^ into the PLM pair with an additional input 
of 2.5 X 10"^ into the ASK sensory neuron pair (known experi¬ 
mentally to promote turning [18]). Simultaneous PLM-I-ASK 
stimulation causes bistability, with relatively long transient 
times r. To the right of each raster plot is the trajectory 
within the Forward-Motion 2D Plane (defined by the trajec¬ 
tory in Panel A, and used for all subsequent projections). 

to construct a full connectomic dynamics model in which 
it was shown to yield reasonable low-dimensional neu¬ 
ral proxies for known behavioral responses (specihcally, 
it was shown that simulating excitation of the tail-touch 
mechanosensory pair PLM creates a two-mode oscilla¬ 
tory limit cycle in the forward motion motorneurons). 
As in m, neural membrane voltage dynamics are gov¬ 
erned by: 

CV) = -G^(U, - E,,u) - /f“'’(V) - /f^”(V) + If"* (1) 

C is the whole-cell membrane capacitance, is the 
membrane leakage conductance and Eceii is the leak¬ 
age potential. The external input current (which we 
change to specify the external stimulus) is given by if'®*, 
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Interaction 

Timescale 

Single-Neuron Membrane Leakage 
Gap Junctions 

Synaptic Connections 

100ms 

10ms 

200ms 


TABLE I. Orders of magnitude for various timescales within 
the system for the parameters chosen. 


while neural interaction via gap junctions and synapses 
is modeled by input currents (gap) and 

(synaptic). Their equations are: 


j 

(2) 

fy- = Y.GtgS,{V.- Ej) 

(3) 


] 


Gap junctions are taken as ohmic resistances connect¬ 
ing each neuron where Gfj is the total conductivity of 
the gap junctions between i and j. Synaptic current is 
proportional to the displacement from reversal potentials 
Ej. Gfj is the maximum total conductivity of synapses 
to i from j, modulated by the synaptic activity variable 
Si, which is governed by: 

Si (Xj-cjy^Vi, K, V)/i)(l S^) 0,fiSi (4) 

where and aj, correspond to the synaptic activity’s 
rise and decay time, and (j) is the sigmoid function 

K, Vth) = 1/(1 + exp{-f3{Vi - Vth)))- 
We keep the parameters values of m- In particular, 
the connectivity parameters and GL are prescribed 
by the full connectome m- The relative significance of 
these specific connectivity values is maintained by not 
htting any of the other global parameters. Instead, these 
parameters are estimated to a reasonable order of mag¬ 
nitude from the literature and assumed equal for each 
neuron. Relevant values to this section are, as taken 
from m- gap junctions and synapses are both given in¬ 
dividual conductances oi g = lOOpS; cell membranes are 
set to a conductance of G° = lOpS; membrane capaci¬ 
tances are set to IpF; and the synaptic rise and decay 
constants are set to = 1 and = 5 s“^. Thus 
all neurons are modeled as identical except for their con¬ 
nectivity and the assignment of them as excitatory or 
inhibitory (where Ej will have one of two values corre¬ 
sponding to these classes). 


B. Model Timescales 

Of particular relevance to this paper are the timescales 
within the system. From the first term in Equation Q, 
we see that the exponential free decay constant of an 
unconnected neuron (i.e. decay through the membrane 
leakage term alone, with 1^°“^ = = 0) would 

be Tfree = GjG'^ = lOOms. Similarly, the time constant 


value given by gap junctions would be Tga-p = Cjg = 
10ms. 

There are also timescales intrinsic to the synaptic 
dynamics. We approximate these by considering the 
dynamics when voltages are held constant, and thus 
4>{vi, K,Vth) = (t>i is constant. Then Equation ( [^ be¬ 
comes: 

Si — ^d)^i (b) 

and thus the synapses will exponentially approach equi¬ 
librium with a time constant of Tsyn = l/(ar/>i + Od)- 
Since = 1 s“^, = 5 s“^, and (pi € (0,1), synapses 
must have exponential time constants in the range Tsyn S 
(166, 200)ms. 

It will be shown that, when the system is in a bistable 
regime, the timescales of transient dynamics within the 
system can be orders of magnitude above any of these 
intrinsic time constants within the system (on the order 
of 10s, for example). 


C. Model Discussion 

The model does not include various extra-synaptic fea¬ 
tures known to drive or regulate responses. For example, 
there is evidence that self-sustained forward locomotion 
in G. elegans is regulated by proprioception within motor 
neurons |21] (compare how our model, lacking this, does 
not sustain oscillation in the absence of explicit external 
input). Computational modeling which includes stretch- 
receptive proprioception shows that such feedback loops 
can control behavioral features such as gait modulation 
between differing environments [22l |23] . The lack of such 
feedback mechanisms and other signaling mechanisms 
(such as various neuromodulators, monoamines and pep¬ 
tides |241 125)1. in combination with the simple neuron 
model and parameter assumptions, mean that specific 
responses to given inputs seen within the model can be 
encoded only within the network’s connectivity. This re¬ 
ductive approach yields information as to how behavioral 
responses could be encoded within the structure of the 
connectome. 

As an example of the model’s ability to generate prox¬ 
ies for behavioral responses encoded within the network, 
it was shown that stimulating the tail-touch mechanosen- 
sory neuron pair PLM, which experimentally leads to for¬ 
ward motion gives rise via a bifurcation to a limit cy¬ 
cle within the forward-motion motorneurons. This limit 
cycle consists of only two modes, in agreement with the 
behavioral observation that the worm’s body shape dur¬ 
ing forward motion is well-described by a similar two¬ 
mode oscillation m- The non-triviality of this agree¬ 
ment was established by showing that simulated abla¬ 
tion studies affected this response in agreement with ex¬ 
perimental ablation studies (e.g., ablation of the AVB 
interneurons destroys the response both experimentally 
and in the model [H])- 
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D. Response of Model to Given Inputs 


Figure shows the response of forward motion mo- 
torneurons to various inputs as a function of time. Panel 
A of the figure shows a raster plot of motorneuron volt¬ 
ages in response to PLM input (through the Text term 
in Equation 0)> for which the two-mode oscillatory re¬ 
sponse can be observed m- The trajectory of these two 
leading modes are plotted as a function of time on the 
right. We use this same low-dimensional space (defined 
as the two forward-motion motorneuron modes which os¬ 
cillate during PLM activation) throughout the paper. In 
other words, we use the same projection for the low¬ 
dimensional trajectories in Panel B and in all further 
figures. 

The response to PLM stimulation alone consists of a 
single possible state (i.e. a limit cycle trajectory), but 
if the model is capable of describing the dynamics in 
terms of long-timescale transitions between states under 
the same input, then we wish to find inputs which allow 
multiple states and transitional dynamics. We find that 
such inputs indeed exist. As an example, we consider the 
response to simultaneous stimulation of the PLM neuron 
pair along with the ASK neuron pair. We choose this 
stimulation since excitation of ASK neurons have been 
shown experimentally to promote turning |18j and their 
ablation greatly increases the duration of periods of for¬ 
ward motion |26j . 

As we show in Figure [ 33 , for this combined input 
there coexist two different attractors, i.e. the system 
is bistable. The two trajectories plotted are in response 
to the same constant input amplitudes into PLM and 
ASK, and differ only by their initial conditions. Note 
that the transients before convergence into the eventual 
fixed point or limit cycle have long timescales (relative 
to the intrinsic timescales of the system as discussed in 
Section IIB). The model therefore does exhibit multista¬ 


bility for this given input, but given the large dimension¬ 
ality of the input space the discovery, identification and 
interpretation of these multistable regimes is not triv¬ 
ial. Since we wish to understand the neural dynamics 
as consisting of long-timescale transitions between dis¬ 
crete attractors, we develop a method for (I) identify¬ 
ing the existence and nature of attractors in response to 
arbitrary inputs, (2) characterizing transient timescales, 
and (3) providing interpretable biophysical meaning to 
calculated trajectories via projection onto a meaningful 
low-dimensional space. 


III. BIFURCATION DIAGRAMS FOR STATE 
IDENTIFICATIONS 

A. Calculating Bifurcation Diagrams 

Motivated by observational studies which describe C. 
elegans behavioral dynamics in terms of low-dimensional 
attractor dynamics |I6j . we wish to understand our sim- 



FIG. 2. The input amplitude into a single neuron at which 
the standard equilibrium becomes unstable varies by neuron 
type. The vertical axis shows the percentage of neurons of a 
given type (sensory neurons, interneurons, or motorneurons) 
for which the standard equilibrium is unstable when said neu¬ 
rons receive the corresponding input amplitude. Note that 
sensory neurons , on average, create a bifurcation at a lower 
input amplitude than do interneurons, which in turn require 
a lower amplitude than motorneurons. 


ulated neural dynamics as consisting of transitions be¬ 
tween discrete attractors. We therefore propose to con¬ 
struct bifurcation diagrams that depict attractors exist¬ 
ing within the system under arbitrary inputs. By fixing 
the direction of the input vector /gajt in Eq. Q and us¬ 
ing its amplitude as our bifurcation parameter, such dia¬ 
grams will show us at a glance the set of states created in 
response to a given input, and provide us with a method 
of identifying induced multistability. 

Figures and 1^ show examples of such bifurcation dia¬ 
grams, in which we plot the furthest distance from 
standard equilibrium (within the 2D Forward-Motion 
Plane) of all attractors present as a function of input 
amplitude. In principle, such diagrams could be calcu¬ 
lated by simply performing a large ensemble of simula¬ 
tions and projecting the results into the 2D plane, but 
such simulations are relatively time-consuming. We can 
take advantage of the fact that, for this model, it is easy 
to compute its Jacobian matrix at any given point for 
any constant input. We therefore use Newton’s method 
when possible, supplementing with simulations to explore 
the input space and find addtitional attractors or when 
Newton’s method does not converge to the desired re¬ 
sult. Full detail on the algorithm used to generate such 
diagrams can be found in the Supplementary Materials. 

We generated these diagrams for all 279 of the single¬ 
neuron inputs into the system. Note that the figures 
from these simulations, as shown in the supplementary 
materials, are done over a much coarser range than those 
in Figures and The purpose of these coarse figures 
is to quickly give an indication of the likely number of 
states for each range of inputs. Thus these diagrams give 
a means of identifying what attractors will exist within 
the system for a broad range of arbitrary inputs, and of 
easily identifying regions of multistability in the input 
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space. 

Generating these diagrams for all possible single inputs 
allows for the qualitative comparison of features within 
each neuron’s bifurcation diagram. Similar features in 
the bifurcation diagrams of neurons may suggest simi¬ 
lar functionalities. As a simple example, in Figure 
we compare the input amplitude at which the standard 
equilibrium first becomes unstable for sensory neurons, 
interneurons and motorneurons. The majority of sensory 
neurons are seen to drive bifurcations in the system at 
lower input levels than for most interneurons, which in 
turn require lower inputs than most motorneurons. In¬ 
tuitively, this suggests that the system is typically more 
sensitive to input into sensory neurons than it is to in¬ 
terneuron or motorneuron inputs. Furthermore, for each 
group of neurons we compute the percentage of single 
neuron inputs which promote limit cycle attractors. We 
find that within our input range, 32.6% of sensory neu¬ 
rons and 26.7% of interneurons give rise to oscillatory 
dynamics, whereas only 8.4% of motorneurons result in 
oscillation when stimulted. This points to the sensitivity 
and particular ability of sensory neurons to drive com¬ 
plex dynamics within the network. Such results serve as 
a demonstration of the ability of these bifurcation dia¬ 
grams to provide meaningful and intuitive information 
about the functionality of neurons and the behavior of 
the system. 


B. A Defining Example: Response to PLM Input 

Figure shows a low-dimensional bifurcation diagram 
for constant PLM input. The figure shows the creation 
of a stable limit cycle in response to input into the neu¬ 
rons PLML/R. By evaluating this bifurcation diagram 
we can identify the regions of interest which have quali¬ 
tatively distinct responses (in this case, the region with 
a lone attractor which is a stable fixed point and the sec¬ 
ond region with a lone stable attractor which is a limit 
cycle after the fixed point attractor becomes unstable). 
For each region we can perform simulations which are 
then projected onto the low-dimensional plane (the PLM 
limit cycle being what defines this plane). Given the 
correspondence of this limit cycle to forward motion, as 
in HU, these low-dimensional trajectories are readily in¬ 
terpretable: the fixed point corresponds to a static worm, 
and the limit cycle corresponds to oscillatory motion of 
the body of the worm. 




FIG. 3. Bifurcation diagram for constant PLM stimulation 
of varying amplitude. Below an input of 1.2 x 10'* the sys¬ 
tem goes to a stable fixed point very close to the standard 
equilibrium, but beyond that input level the system goes to a 
stable limit cycle (where the plotted point gives the furthest 
distance from standard equilibrium on the limit cycle). The 
diagram shows the two qualitatively distinct regions of inter¬ 
est for PLM inputs: the low input level in which the system 
remains at a fixed point, and the higher input level beyond 
which the system enters into a limit cycle (which in this case 
can be considered to serve as a proxy for forward motion m)- 


the input into the ASK pair. Figurej^shows the resulting 
bifurcation diagram. At inputs below 1.5 x 10'^, the limit 
cycle remains relatively undisturbed. At greater inputs, 
however, a series of bifurcations occur such that there 
is a sudden jump in the distance of the limit cycle, and 
at about 1.7 x 10"*^ the system becomes bistable with the 
addition of a new fixed point. Thus we are able to imme¬ 
diately identify from this figure multistability within the 
system, which we may then go on to investigate further. 
Specifically, we are interested in the further investigation 
of transient timescales of the system. 


C. Characterizing Bistable Dynamics 


Of greater interest are responses to compound activa¬ 
tions; that is, more complicated inputs leading to more 
complicated responses. We consider as an example the 
dual input into the P LM and ASK neuron pairs as dis¬ 
cussed in Section IID We keep a constant input of 2 x 10^ 


into the PLM pair and use as our bifurcation parameter 


IV. LONG TIMESCALES AND INTERSTATE 
TRANSIENTS 


In Figure we investigate spatial and temporal as¬ 
pects of the convergence onto one of the two bistable 
states. An ensemble of 200 simulations (with random ini¬ 
tial conditions in the neighborhood of the standard equi- 
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FIG. 4. Bifurcation Diagram for varying amplitude of input 
into the ASK pair. Input into PLM is fixed at 2 x 10"^. Note 
that as input into ASK increases, the forward-motion limit 
cycle remains relatively undisturbed until it reaches about 
1.5 X 10"^, after which the distance jumps and a fixed point 
becomes stable, giving rise to a bistability within the system. 


librium) were performed for each ASK input level. From 
those, the trials converging to the fixed point solution 
were taken and the convergence time t was calculated 
by calculating, for each fixed point trial, the time after 
which all points of the trajectory are within a distance 
e of the final value (using here e = 0.004). The average 
and standard deviation of these convergence times are 
shown in the top right of Figure Convergence times 
for the limit cycle solutions are qualitatively similar when 
comparing trajectories such as those in the upper-left of 
the figure. Note that these convergence times are con¬ 
siderably longer than other timescales within the system 
(comparing, for example, the model’s free neuron decay 
constant of 10ms [13], or other trajectory timescales such 
as the limit cycle period, which remains approximately 
two seconds regardless of ASK input). 

Shown also are the basins of attraction for trajecto¬ 
ries starting on the low-dimensional plane, on a grid of 
initial conditions centered at the standard equilibrium 
point (which we choose as our origin). The size of the 
grid is chosen to be within a small neighborhood of zero 
(within the range (—4,4) x 10“®) since we hnd that tra¬ 
jectories initiated farther away are hrst attracted towards 
the zero point before being rerouted to the fixed point or 
limit cycle attractors. We observe consistency between 
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FIG. 5. Spatial and temporal properties of convergence 
for PLM-FASK input (i.e. the bistable region of Figure]^. 
The upper-right plot shows fixed point convergence times as 
a function of input amplitude (from 200 trials at each point). 
Note the relatively long transient timescale. The second row 
shows the spatial basins of attraction for different inputs. 
Each grid covers a small region around the standard equi¬ 
librium, plotting on (—4,4) x 10“® for both modes. At an 
ASK input of 1.6 x 10"^ all initial conditions converge to a 
limit cycle, but initial conditions on the plane are split be¬ 
tween the limit cycle and fixed points at higher inputs such 
as 2.4 X 10"*. 


the structure of these basins of attraction and the ob¬ 
served convergence times. In addition, the basin of at¬ 
traction plots indicate distinct regions in which initial 
conditions are more probable to be attracted to the fixed 
point, and thus shows which portions of forward move¬ 
ment (i.e. which segment of the PLM-driven limit cycle 
trajectory) are more prone to ASK-driven transitions. 


V. CONCLUSION 

We explored the input space of a C. elegans neural 
dynamic model which incorporates its fully-resolved con- 
nectome and demonstrated that various multistabilities 
arise in response to inputs. Using a low-dimensional pro¬ 
jection space based upon forward motion, we are able 
to systematically explore responses to complex inputs 
and understand them in a framework of low-dimensional 
attractor dynamics. In our study, the bifurcation dia¬ 
gram is constructed by using the constant-in-time input 
as our bifurcation parameter. We show that such di¬ 
agrams are capable of revealing and mapping multiple 
attractors within the system by using a low-dimensional 
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projection space which guides the search for attractors, 
identifying their stability and their effect upon forward 
movement. Furthermore, the low-dimensional projec¬ 
tion helps in the interpretation of the dynamics upon 
the discovered attractors, especially the dynamics asso¬ 
ciated with multistability. We characterize such multi¬ 
stable dynamics, noting specifically that when the sys¬ 
tem enters into a multistable regime, transient timescales 
within the system can be very long relative to intrinsic 
neural timescales (comparing, for example, the three or¬ 
ders of magnitude between the ~ 100ms neural timescales 
in Section IIB to the ^lOs transient lengths in Figure]^. 


The fact that multistability within the connectomic dy¬ 
namical system is capable of generating such long tran¬ 
sient timescales has critical implications. These longer 
timescales, on the order of many seconds, are on a simi¬ 
lar order to many behavioral timescales such as forward 
crawling survival time m- This suggests that various 
behavioral responses could be associated not with the 
attractor itself, but rather with the transient leading to 
that attractor. This is consistent with theoretical con¬ 
structions and experimental observations of transient or¬ 
bits between attractors m- Importantly, this view¬ 
point is supported independently and in a completely dif¬ 
ferent theoretical framework by direct connectomic sim¬ 
ulations from biophysically appropriate neuron dynamics 
within the worm, i.e. the multistability of attractors and 
long-time transients are not engineered in the model to 
fit the data and observations, rather they naturally arise 
from the dynamics associated with the connectome. 


This study suggests that neural computations can con¬ 


sist of both dynamics on attractors (as in our PLM-driven 
limit cycle) and of long-timescale transients between mul¬ 
tiple attractors which may arise in the system (as we show 
in the long-timescale transients between the multistable 
states from PLM-I-ASK input). We have demonstrated 
that both dynamical features can arise by applying sim¬ 
ple, identical neuron models onto the C. elegans connec¬ 
tome data, suggesting that these responses are encoded 
within the connectome itself. 

More broadly, many networked dynamical systems 
across the engineering, physical and biological sciences 
may also be dominated by patterns of activity and long¬ 
time transients induced by the structure of the network 
architecture. Understanding the basic principles of such 
behaviors is critical for optimizing performance and con¬ 
trolling deleterious effects. One can easily imagine sce¬ 
narios in which suppressing transients would be impor¬ 
tant, such as the observed power-grid swing instabil¬ 
ity m- The analysis above may be able to help under¬ 
stand how the network architecture encodes such dele¬ 
terious patterns of activity when combined with rele¬ 
vant dynamics. In contrast, one might desire to gen¬ 
erate a network architecture to induce a transient that 
is beneficial for some purpose relative to an application 
(for instance, a crawling motion in the case of the C. 
elegans). Understanding how the network connectivity 
graph drives such activity would be critical for inducing 
such beneficial patterns of activity, perhaps even sug¬ 
gesting network control protocols for achieving desired 
results. The theoretical framework presented here high¬ 
lights the rich and complex dynamics that emerge with 
networked architectures. 
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